


GENERALIZED SELECTIVE MODAL ANALYSIS 

JULIAN BARQUiN * 

Abstract. A new approach which generalizes the Selective Modal Analyis (SMA) and algorithms 
based upon it for solving the generalized eigenvalue problem is described. This approach allows for 
the systematic consideration of physical properties of the system under study. Two small application 
cases demonstrate the capabilities of the proposed approach. 
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1. Review of SMA. Selective Modal Analysis (SMA) is a physically motivated 
framework for understanding, simplifying and analyzing complex linear time invariant 
models of dynamic systems. SMA can focus on selected portions of the structure and 
behaviour of the system 0, ||] . 

In many physical systems, it can be readily recognized that some set of modes 
is asociated to a certain set of variables. For instance: in electrical power systems, 
the electromechanical oscillations are associated to the machines rotors's angle and 
speed. Very often, this association is used, explictly or implicitly, to simplify the 
mathematical models of the system under study. SMA aims to exploit this relationship 
in a systematic and rigorous way. 

So, let us assume that it is desired to analyze the dynamic system: 
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The SMA approach is to classify the x components in relevant (r) and less-relevant 
(z) components. So, possibly after a trivial reordering, it can be written: 



(1.2) 



Therefore, equation ( n^ ) can be written as: 



(1.3) 



A, 



A, 



Let us assume that there is interest in computing and eigenvalue A, and its left 
and right eigenvectors v and w. 



(1.4) 
(1.5) 



Aw — Av 
u^X = w^A 



The eigenvectors can be partitioned analogously as the states x: 



(1.6) 
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It is easy to check that it must be fulfilled that: 



(1.7) 
(1.8) 



XVj. — (j4,.j. + Arz\X ~ Azz) Azr) Vr 
wJX = wj (Arr + Arz{X " Azz^^Azr) 



So, the interesting eigenvalue A is in the spectrum of the matrix Arr + Arz{X ~ 
Azz)~^Azr- On the other hand, if this mode is strongly correlated with the relevant 
variables r, it should be expected that the spectrum of Arr contains an eigenvalue 
quite similar to A and, therefore, that the matrix Arz{X~ Azz)^^ Azr perturbs sligthly 
the desired mode. That suggests the following algorithm: 



Algorithm 1 






Input : Arr ,Arz, Azr ,Azz- 






Output: X,Vr,Wr- 






1. 


Perform the eigenanalysis of Arr, 






2. 


Select the interesting mode °A, "w^, ° 


Wr, 




3. 


for j = 1, 2, 3, . . . until convergence. 








3.1. Compute H{^-^X) = Arzi^^^X 


-Azz)- 


'Azr. 




3.2. Perform the eigenanalysis of Arr + H{^ 


-'X), 




3.3. Select the interesting mode ■'A, 


'Vr,^Wr, 




4. 


end 







The convergence properties of this algorithm have been studied in |3| . The conver- 
gence is controlled by the eigenvectors. Especifically, the algorithm locally converges 
if and only if it is fulffiled: 



(1.9) 



P 



wjvr 



W^Vz 



> 1 



The number p is called the participation ratio. Notice that, for any eigenvalue, 
this ratio only depends in the way that the variables have been partitioned in relevant 
and less relevant. 

The former algorithm can be generalized to search for several eigenvalues. The 
interesting eigenvalues shall be collected in a diagonal matrix A, and the eigenvectors 
in matrices V and W: 



Ai 
A2 



(1.10) A= . ." . . V ^ [VIV2 . ■ ■ Vr,] W ^[wiW2...Wn] 

... A„ 
The matrices V and W can be also partitioned in relevant and less relevant parts: 



(1.11) 



V = 



Vr 
Vz 



w = 



Wr 

Wz 
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So, the following algorithm can be proposed: 



Algorithm 2 


Input : Arr ,Ar^,A„,A^^. 


Output: K,Vr,Wr. 


1. Perform the eigenanalysis oi A^ 


2. Select the interesting modes ^Aj^K-i^Wr, 


3. for j = 1, 2, 3, . . . until convergence, 


3.1. Compute a matrix ^ M which fulfills 


np-^Vr= [H{^-^Xiy-^Vir,Hi^-^X2y-^V2r, ■ ■ ■ Hi^-^X^y-^Vnr] 


where H{X) = Ar,{X - A,,)-^A„, 


3.2. Perform the eigenanalysis of A^ + ■' M, 


3.3. Select the interesting modes ^ K^Vr^Wr, 


4. end 



The convergence conditions of algorithm 2 are also studied in H , although they 
are considerably more involved than those of algorithm 1. However, the computa- 
tional experience shows that the convergence is good if the participation ratios of the 
interesting modes are high. 

From the point of view of the computational effort, the most demanding task of 
both algorithms is the step 3.1, which requires to solve linear system involving the 
matrix X — Azz- Most SMA applications have been developped for the study of electric 
power systems |^, where special techniques based upon the peculiar characteristics 
of these systemas have been used to perform efficiently this task. 

In addition to algorithms 1 and 2, there are a number of related ones which 
considerably improve their speed and robustness |^, |6|. However, these algorithms 
are outside of the scope of this paper. 

2. Generalized SMA. Although SMA has been succesfuUy used in applica- 
tiosn, it has some drawbacks: 

1. Very often, the computation of the desired modes begin with the study of 
a simplified model, in order to gain an insight on the interesting eigenstructure. Al- 
though SMA allows to use subsequently the information on the relevant variables, it 
does not make use of all the information obtained with the simplified model. 

2. There are some problems, specially in continous media, where it is known the 
overall shape of the desired modes, but it is impossible to adscribe them to a small 
number of system variables. 

Both kind of problems are treated in this paper examples. The aim of this section 
is to generalize the classical SMA theory in order to deal with these problems. 
So, the problem to solve is the eigenvalue problem: 



(2.1) XEv = Av 

(2.2) w'^EX^w'^A 

A and E are mxm real matrices. The matrix _E is a symmetric, possibly singular, 
projection matrix: 
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(2.3) E^ =E = E^ = E'^ 

The superscript T denotes the transpose and f the hermitian conjugate. It is 
assumed that the right eigenvector v approximately hes in the subspace spanned by 
{ei, 62, ... , e„}. UsuaUy, n <^ m. Then, it is defined the matrix 



(2.4) £:=[ei,e2,...,e„] 

In similar way, the left eigenvector w is assumed to yield, approximately, in the 
subspace spanned by {/i, /2, . . . , /«}• So, it is defined the matrix 



(2.5) ^=[/i,/2,...,/„] 

Besides, the e^ and fj basis are normalized in order to fulfill the equation: 



(2.6) Pes ^ In 

/„ is the nxn identity matrix. This equation can be enforced so long as no vector 
generated by the basis e^ or the basis fj is included in the kernel of E. This condition 
shall be assumed in the sequel. 

Then, the eigenvectors v nd w can be written as: 



(2.7) v = £a + z 

(2.8) T'^Ez == 

(2.9) w = TI3 + y 

(2.10) S'^Ey^Q 

It is easy to show that the above decomposition exists and is unique. Then, after 
the algebraic manipulations shown in the appendix A, it is found that: 

(2.11) \a^ Arra + H{\)a 
where 



-1 , 



(2.12) Arr = T'^AE 

(2.13) H(\) =T'^AV{\E~A+\A,(3\^Y'VAE, 

The matrices Q and V are idempotent matrices defined by 

(2.14) V = /,„ - EET^E = Ira - Q 
and [Q, ^] I is the anti-commutator: 
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(2.15) 



UAl 



QA + AQ 



These formulae are the basic ones in selective modal analysis, and can be consid- 
ered as a generalization of "classical" SMA, as shown in the appendix B. It is also 
easy to check that (3 is the left eigenvector of Arr + H{\). Notice that the dimension 
of the matrices A^ and H{X) is n <^ m. 

It is also noteworhy that 



(2.16) 
(2.17) 



z = V{XE-A+[A,Q]^} ^VASc 
yt ^ p^j^^AV {XE-A+ [A, Q]+Y' V 



, -1 , 



3. Algorithms based on H{\) computation. The aim of this section is to 
apply the results of the former section in order to obtain workable algorithms, and 
their convergence conditions. In order to simplify the notation, let us denote by 
N{\, Q) the matrix 



(3.1) 
so that 



N{\,Q.)={\E-A+[A,Q] + y 



(3.2) H{\)=T^AV{^-^\E-A+[A,Q]^] ^ VA£ ^ T^^ AVN{\Q)VAS 

3.1. Linear algorithm. Specifically, let us consider the following generalization 
of algorithm 1 : 



Algorithm 3 
Input: E,A,£,T. 
Output: A,a,/3. 

1. Form Arr — ^'^ A£, and perform the eigenanalysis of Arr, 

2. Select the interesting mode ''A,''a, "/?, 

3. for j = 1, 2, 3, . . . until convergence, 

3.1. Compute H{^-^X) = T'^AVN{^-'^\ Q)VA£. 

3.2. Perform the eigenanalysis of Arr + H{^^^X), 

3.3. Select the interesting mode ^XjaJ (3, 

4. end 



The following theorem states the conditions for the local convergence of Algorithm 
3: 

Theorem 3.1. Given an eigenvalue A of the pair {E,A) with associated right 
and left eigenvectors v and w, there is 5 > such that if ||°A — Aj| < S, algorithm 3 
converges to the eigenvalue A whenever 



(3.3) 



\P\ 



w'^ Qv 
w'fiE- Q)v 



>1 
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Furtthermore, the error ^e — ^ X — X fulfills: 



(3.4) 



oJ-i. 



o{- 



i-K 



Proof. The proof is given in the appendix C. D 

It is eas y to check that p just defined is, in the "classical" case, the same p defined 
in equation (|l.9D. 



3.2. Superlinear algorithm. If, in algorithm 3, the right eigenvector v (respec- 
tively the left eigenvector w) is contained in the span of £ [T) , then z = {y = 0) 
and p ~* oo. So, the algorithm could be speeded up if the matrices £ and J-" are 
pdated in order that their span contains the last approximation to v and w. ^v and 
^w. Therefore, the following algorithm is proposed: 



Algorithm 4 








Input: 


E,A,°8,'>J^. 








Output: X,a,P,v,w. 








1. Form ^Ar-r = °JpU °£, and perform 


the 


eigenanalysis of Arr, 


2. Select the interesting mode °A, °a, ° 


/?, 






3. for 


7 = 1,2,3,... until convergence, 








3.1. 


Compute Hi^'^X) ^ 










j-ijrt^j-i7?Ar(i-iA,J-iQ) i- 


^VA 


J"-i£. 




3.2. 


Perform the eigenanalysis of ^~ 


1 A 


+ H{i- 


^A), 


3.3. 


Select the interesting mode ^A, 


^a,i 


/?, 




3.4. 


Compute ^z = ^'^VN{^~^X,^~ 


'Q) 


i-i-pA^- 


-^£^a 




and Jyt = ^ f3l3-^Tl A3-^VN{^ 


''X, 


;-ig).- 


-'V, 


3.5. 


Compute ^v = ^^^E^a + ^ z and ■'«;''' 


= J'/^tJ 


-i.Ft+^yt, 


3.6. 


Update ^£, ^ T in such a way that 








^v e span(-'£),-'u' e span(^J^), 








3.7 


Form^r^ ^^T'^A^e, 








4. end 











A particular case of algorithm 4 is when the matrices ^£ and ^T are vectors. 
Then, these matrices are esentially the estimated eigenvectors. The local convergence 
properties of the algorithm are, in this case, given by the following theorem: 

Theorem 3.2. Given an eigenvalue X of the pair {E,A) with associated right 
and left eigenvectors v and w, if \\N{X, Evw"^ E)\\ < cxi, there is a neighborhood of 
X,v,w, such that i/°A,'^w,°it; belong to it, algorithm 4 converges. Furthermore, it is 
fullfilled that asymptotically there is a constant K .such that 



(3.5) 



< KP^\ 



|l + \/2 



Proof The proof is given in the appendix ^. D 

In the general case, whenever ^£ and ^J- are not vectors, it is expected that 
algorithm 4 converges at least so fast. This is because the relevant subspace is bigger, 
so that the approximation to the eigenstructure can not be worst. 
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3.3. £ and T selection. The computation of H{\) requieres to make a selection 
of the matrices £ and J- . The foUowing theorem can be used for this task: 

Theorem 3.3. 

The matrix A„-\-H{X) is invariant under the transformations £ ^ £ + (I„i — E)C 
or J-' ^ J-' + (/,„ — E)Ai, where C or Ai are arbitrary matrices of the same dimension 
than £ or T . 

Proof. 

The proof is provided in appendix O 



D 



In many cases the E matrix can be written as: 



(3.6) 



E 



Ir 




Ir is a r-dimensional identity matrix. Therefore, the vectors and matrices can be 
partitioned in dynamic and static parts. Specifically, 



(3.7) 
(3.8) 



£ 



£d 
£s 



The invariance of the matrix A„ + ^(A) under the considered transformations 
means that the value of the static components £s and Tg is irrelevant in order to 
compute this matrix. 

So, referring to algorithm 4, there are at least two possibilities: 

1. To keep the whole eigenvector in step 3.4: update -^S, ^T in such a way 
that ^v G span{^£),^w G span(''JF). Then, as the algorithm converges, the matrix 
H{X) -> 0. 

2. To update the matrix as above, but the static components, which are zeroed 
{£s = and Ts = 0). As the matrix converges, the matrix H{X) converges to a 
non-zero value. 

The second possibility can be useful in order to minimize the numerical effort. 

4. H{X) computation. From the computational point of view, the most de- 
manding task of the algorithm is the computation of the matrix H{X) or of the matrix 
Arr + H{X). The purporse of this section is to propose algorithms to deal efficiently 
with this task. 

4.1. H{X) computation using the Shermann- Morrison lemma. H{X) can 
be written as: 



(4.1) H{X)=^J'^Ar{XE-A+[A,Q]^} rA£ 

The basic problem is related to the matrix 



(4.2) 



N{X)^{XE-A+[A,Q]^} ' 
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A problem is that, generally, the matrix [A, Q], is not sparse. However, it is 
possible to obtain an expression of N{X), which allows computations by using only 
sparse matrices, by means of the Shermann-Morrison lemma. So, 



(4.3) 

(4.4) 

(4.5) 
(4.6) 

(4.7) 
(4.8) 



= hE-A+ [A,E£T^E]\ 

= {XE-A+ (AES) [T^E) + {E£) {pEA)y^ 
= {XE-A+ {AE£) {Pe) + (ES) {T^EA)y^ 
= {XE - A + Tjcj)^ - r)(t)^ + 

{AE£) {Pe) + {E£) {T^EA)Y^ 



XE-A + rj4)'^-[ f] -AES 



-ES 



T^E 
T'^EA 



T] and (f) are two sparse vectors which make sure that the matrix XE — A + rjcjr 
is regular even if A is an eigenvalue of the pair (A, E) . This is going to happen when 
the SMA algorithm converges. Now, from the Shermann-Morrison lemma: 

N{X) = {XE-A + r^(t)'^y^ + 

{A^-A + 770t}"^ [ r/ -AES -ES ] 



I: 



(4.9) 



2n+l 



T^E 
T^EA 



T'^E 
T'^EA 



{XE- A^ri(l>^) ^ [ ^ 



-AES -ES 



{XE- A^^<^'^\ 



The number 2n + 1 is usually small. So, it is only required to know the LU 
factorization of a filled matrix of small (2n + 1) dimension and of the sparse matrix 
XE- A + ri^'^. 

4.2. Composite models. In power systems analysis, the system to analyze is 
a set of dynamical subsystems connected through a static relationship. Specifically, 
there are I subsytems 



(4.10) 
(4.11) 



EkXMk — AkXMk + BkXjk 
XOk = CkXMk + DkXjk 



1,- 



J 



The variables XMk are the state variables of the fc-th subsytem, xik are the input 
variables and xok the output variables. It is assumed that the number of output and 
input variables of each system is equal. It is also assumed that the matrices E^ are 
sy metric projection real matrices: 



(4.12) 



Eh 



El 



E. 



El 
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Let US define the vectors 



(4.13) 



Xm 



xmi 




Xll 




xoi 


XM2 


XI = 


XI2 


xo = 


xoi 


XMI 




Xll 




xoi 



In addition to these equations, there is also a static interconnection: 



(4.14) 



Jll Jl2 
J2I J22 



XI 
XA 



Xo 




XA is a set of additional algebraic variables. 

For this kind of systems, it is convenient to consider the following £ and T 
matrices: 



(4.15) 



which mimics the x structure. Furthermore, £m and Tm are defined on a subsytem 
basis: 



£m 




T = 


^M 


















£mi 


.. 


. 


(4.16) 


£m = 





£m2 ■ ■ 


. 









.. 


• £mi 






^Ml 





.. 


(4.17) 


^M = 





^M2 ■ 


.. 


mi 






. 1 • l^ 








diag{£Mi,£M2, ■ ■ ■ , £mi) 



diag(JFA/i , ^M2 , ■ ■ ■ , ^Mi ) 



Then, as proved in the appendix O, the matrix H{\) can be computed as 



H{\) ^ Ha{\) + [{Br + Hb{X)) 0] 



Jii-{D + Hd{\)) J12 
J21 J22 



n -1 



(4.18) 

where all the matrices are computed in a subsytem basis: 



Cr+Hc{\) 





(4.19) 
(4.20) 
(4.21) 
(4.22) 



Ar = dl&g{Arl ...Arl) 

Br = diag(Bri . ■ . Bri) 
Cr = diag(Cri . ..Cri) 
D = dia.g{Di . . . Di) 



10 J. BARQUIN 

(4.23) HA = diagiHAi...HAi) 

(4.24) Hb = dis.g{HBi . . . Hbi) 

(4.25) He = dmg{Hci . . . Hci) 

(4.26) Hd = diagiHoi ■ ■ ■ Hdi) 

5. Direct algorithms. As said above, the most difficult task in order to apply 
SMA is the computation of the matrix H{X). However, the proposed algorithms can 
be formulated without needing to compute this matrix. The purporse of this section 
is to explain the way of doing it. 

5.1. Single eigenvalue algorithms. The basic SMA formulae are: 



(5.1) ^X^a = {Arr + H{^-^\)) ^a 

(5.2) ^'/gtJA = ^'/Jt {Arr + H{^-^X)) 

(5.3) ^z = V {^-^XE -A+[A, Q] + y^ VAS^a 

(5.4) Jyt = f3^j^fAV {'-^XE -A+[A, Q] + }"' V 

as shown in prvious sections. Let us also, as above, define the vectors 

(5.5) h> = £^a + ^z 

(5.6) ^w = T^/3 + ^y 

Then, the following theorem can be stated: 

Theorem 5.1. Let us assume that equations (5.1-5.4) are fullfiUed. Then, the 
vectors ^v and ^w fullfill the equation 



(5.7) [A-^-'^XiE-Q)]^v^^XQ^v 

(5.8) ^u>t [A - J'-^A {E ~ Q)] = ^w^Q^A 

Proof. The proof is provided in Appendix ^ D 



Let us write equation (5/7) in the following equivalent form 



(5.9) (A - ^-^XE) H = (J A - ^-^A) Q^v 
As T'^ E^ z — 0, it is obtained 

(5.10) Q^v = E£^a 
Therefore 

(5.11) (A - ^-^A^;) ^v = (^A - ^-^A) EE'a 
Let us define 
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(5.12) 



So 



IV ^ {A-^-^\E) ^ EE 



(5.13) 

Analogously, 



3v^^V^a{^\~'-^X) 



(5.14) 
(5.15) 



'W^ =T^E{A~'-^XE) ' 



Now, let us consider the matrix ^W^ (A — ^ ^A£') -'V. 



(5.16) 



^W^ {A - ^-^XE) ^V^a = T^E^v pA - ^^^A)" 

= {^x-^-^xy^^a 



Therefore, ■'a is a right eigenvector of this matrix with associated eigenvalue 
(^A — ■'"^A) . It can be easily shown that ^ j3 is the left eigenvector. Note that, 
because of ( ^.14 ), this matrix can be also writen as T^E^V . 

Equations ( j.lSJ ), (5.14) y (5.16) provide an alternative way of performing the 
SMA iteration. Specifically, algorithm 3 can be also writen as: 



Algorithm 5 






Input: E,A,8,T. 






Output: A,a,/3. 






1. Form Arr = T'^ AS, and perform the eigenanalysis of A^, 

2. Select the interesting mode ^X,^a,^f3, 

3. for j = 1, 2, 3, . . . until convergence. 


3.1. Compute ^Y = 

3.2. Compute m^ 

3.3. Perform the ei 


= {A-i-^XE) E£, 
^J'^E{A-3-^XEy\ 
genanalysis of ^7W = T^iE^V 




3.4. Select the interesting mode ■'a,-'/?, (^A — ■'^^ 

3.5. Update ^A, 
4. end 


A)-\ 



As this algorithm is esentially the same one that Algorithm 3, its convergence 
conditions are the same ones. Note also that step 3.2 could be omitted. Analogously, 
Algorithm 4 can be written as: 
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Algorithm 6 








Input: 


^,A,"£,"^. 








Output: X,a,f3,v,w. 








1. Form °Arr = °J^U 


^£, and perform the 


eigenanalysis of ^Arr, 


2. Select the interestin 


gmodeOA,*'a,°/3, 






3. for . 


7 = 1,2,3,... until convergence, 






3.1. 


Compute W = 


{A-i-^\E) El 


-^£, 




3.2. 


Compute ^W^t 


^i-^T'^EiA-i-^ 


\E) 


-1 

1 


3.3. 


Perform the eigenanalysis oi ^M = 


j-i 


J'^EW, 


3.4. 


Select the interesting mode ^aj (3, 


(n- 


-i-'xy\ 


3.5. 


Update ^A, 








3.6. 


Update ^£, ^T 
^Yespan(J£), 


in such a way that 
JW^espan(J'j^), 






4. end 











Step 3.6 warrants that ^v € span(^£),-'u; g span(^jr). The simplest way to 
achieve that is to set i£ = ^V and ■'JF — ^W but, possibly, a normalization constant. 
As previsously as Algorithm 5, convergence results related to Algorithm 4 can be 
directly applied to Algorithm 6. 

5.2. Several eigenvalues algorithms. In this section algorithm 2 is written 
in the direct formulation, and a new one which includes relevant subspaces updated 
is also proposed. However, no convergence results are provided. 

Let us firstly define the natrix 

(5.17) ^V ^[A-i-^X{E-Q)Y^ EE 
It is clear, from (^ and (FTo|) : 

(5.18) H = ^XiV'a 
Analogously, it is possible to write: 

(5.19) 'W^ ^T^E[A-'-^X[E-Q)Y^ 

(5.20) '^w^ = ^pm^^X 
Let is consider the matrix 

(5.21) W"^ = W^ [A - ^-^X [E ~ Q)] "^ J'F = Pe^V = m'^EE 
It is fuUfiUed, because of (ISASi) : 



(5.22) 
(5.23) 
(5.24) 



= ^X-^T^E^v 
= ^X-^'a 
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Analogously, 



(5.25) 



J"/3tJAA-i = J/jtJA-i 



Therefore, the matrix ^ M contains the sought eigenvalue in its spectrum. On the 
other hand, it can be shown (see appendix |) that 



(5.26) 



^n = ^M-^ + ^-^\in 



where /„ is the identity matrix. Besides, from (5.1) and (5.2), 



(5.27) 
(5.28) 






In fact, a stronger result can be obtained: 

Theorem 5.2. H{X) = M - Arr 

Proof. The proof is provided in the appendix [j|. D 

Therefore, it is proposed the following generalization of algorithm 2: 



Algorithm 7 






Input: E,A,£,T. 






Output: \k,ak,f3k- 






1. Form Arr — J^^'AS, and pei 


form the eigenanalysis of Arr, 


2. Select the interesting modes °Afe,°afe,°/3fc, 




3. for j = 1, 2, 3, . . . until convergence. 




3.1 For each eigenvalue k = 


--1,...,K, 




3.1.1. Compute ^Vk = 


{A-^-^\kEy 


'^ E£, 


3.1.2. Compute ^Mk = 


= T^E^Vk, 




3.1.3. Compute Wfc = 


^M^^+'-^Xk 


In, 


3.1.4 Compute ^hk = 


(Wfc - Arr) '- 


'ak. 


3.2 Compute a matrix ^M which fuUfills 




ni['-^ai,...,^-'^aK 


= '^h^,...,{hK\ 


3.3 Perform the eigenanaly 


sis of Arr +^M, 


3.4 Select the interesting modes ■'Afe, ^ak, 


'Pk, 


4. end 







This algorithm is nothing else that the direct version of algorithm 2, and it reduces 
to it if the S and T matrices are chosen as shown in appendix H. A superlinear version 
of this algorithm, by updating the £ and J-" matrices, is also proposed: 
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Algorithm 8 
Input: E,A,^E,^T. 
Output: \k,ak,Pk- 



Form ^Arr = '''J-'^A'^E, and perform the eigenanalysis of '^A^. 
Select the interesting modes °Ai;,°afe,°/3fe, 
for j — 1,2,3, . . . until convergence, 
3.1 For each eigenvalue k — 1, . . . ,K, 

3.1.1. Compute Wk = {A - ^^^XkE) 

3.1.2. Compute ^Wl ^ ^-^T^E {A - 

3.1.3. Compute ^ Mk ^ T^ EWk, 



-^XkEY 



3.1.4. Compute Wfe = ^Mu^ +^-^\kI, 



3.1.5 Compute ^hk = (Wfe - •''"Mr. 
3.2 Compute a matrix ^M which fuUfills 



J-1 



afc, 



JM p-i 



Qfl, 



i-1 



"K 



p/ii, 



3.3 Perform the eigenanalysis of ^ ^ A, 



3.4 Select the interesting modes ^Xk, ■'oik, -' Pk, 

3.5 For each eigenvalue k = 1, . . . ,K, compute 



'Vk - 



'Vk^ak, and 



'pVw! 



fe' 



3.6 Update ^£ and ^JF in such a way that 



'Vl,. 



^vk] e span(-'£), 



'wi,. . . ,-'wk] e span(-'.F), 

3.7 Compute ^ Arr = 'T^A^E, 

3.8 Compute ^ ak and ^ Pk by imposing 



^E^ak = ■'wfe and ^J^-'/^fc = ■'w^, 



4. end 



Step 3.5 requires to pair each eigenvalue obtained of the eigeanalysis of Arr +^ M 
with the eigenvalues of the previous iterations. A way to do that is to pair trying to 
maximize the scalar products ^dl^^ak2- 

6. Numerical tests. The aim of this section is to apply Generalized SMA to 
two very different problems: the computation of the natural frequencies of a plate 
and the computation of the electromechanical modes of an electric power system. 
The code was developed in MATLAB language. 

6.1. Natural frequencies of a cross-shaped plate.. It is intended to com- 
pute a natural frequency in a cross-shaped plate with unequal arms. Mathematically, 
the problem to solve is 



(6.1) 



- Alp = Lu^ip tl){x, y) e il 



where fl is the cross-shaped dominion, lo the natural frequency and -ip the sought 
mode. In addition, the following boundary conditions must be fulfilled: 



(6.2) 



'ipix,y) =0 y{x,y) e dn 



GENERALIZED SELECTIVE MODAL ANALYSIS 



15 





Figure 6.1. Initial estimation and computed mode 



The previous partial differential equation is approximated by a finite differences 
sdieme. In order to apply algorithms 5, 6, 7 and 8, it is decided that: 

1. As it is intended to solve a symmetrical problem, matrices £ and J-' are taken 
to be real and equal. Furthermore, the hermitian operator f can de substituted by 
the transpose operator T, and the algorithms can be programmed in real, instead of 
complex, algebra. 

2. In algorithms 5 and 7, matrix £ is a vector which approximates the sought 
mode. In algorithms 6 and 8, matrix £ is initially that same vector, and it is updated 
in each iteration to the last mode estimation ^v. 

It was decided to compute the mode corresponding to the upper arm of the cross 
oscillating against the lower one. The initial mode estimation for algorithms 5 and 6 
is shown in figure (6.1), as well as the computed mode. Notice that the algorithms 
converge to the mode whose shape is closest to the initial assumption, being any initial 
assumption of the value of the sought eigenvalue largely irrelevant. 



The eigenvalue is A = ijj^ = 0.1157. Figure 6.2 shows the evolution of the absolute 
value of the error ■'A — A for both algorithms 5 and 6. 

To apply algorithms 7 and 8 an additional mode w^ = 0.1243, co rresp onding to 
the right arm oscillating against the left arm, was computed. Figure |6.3| shows the 
evolution of the absolute value of the error (algorithm 7 in solid lines and algorihtm 
8 in dotted lines). 

6.2. Electromechanical modes of an electric power system. The electric 
power system represented in figure 3.4 (a simplified model of the New England electric 
power system) was analyzed by using SMA. The circles reprents electric generators 
and the lines the electric transmission lines. The electric generators are modelled 
by 9th to 11th order linear systems, whilst the electric network is modelled as an 
algebraic constraint. Therefore, the system is a composite system, as explained in the 
previous section. 

The generators input variables xjk sue the axis and quadrature terminal current 
(two variables per generator), and the output variables xok the axis and quadrature 
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15 20 25 

# of iterations 



Figure 6.2. uj convergence. Algorithms 5 and 6. 

terminal voltage (two variables per generator). The algebraic variables xa are net- 
work voltages. The state variables XMk are mechanical, electromagnetic and control 
generator variables. 

The following physical information is known a priori: 

1. The most troublesome modes are those called electromechanical modes. These 
modes arc related to the generators angle 5k and speed oJk- It is fuUfilled that 
5k = l20TTLUk. The relationship of uJk with the rest of the variables is much more 
complex. However, in a very rough approximation, the subsystem {5k,ujk) can be 
considered a damped pendulum. 

2. It is known that the frequencies of the electromechanical eigenvalues are in 
the order of 1 Hz. 

3. The electromechanical modes can be understood as oscillations of one gen- 
erator or group of generators against other generator or group of generators. 

Fom fact 1, it follows that the 5k and ujk right eigencvector components fulffiU: 



(6.3) 



^v{5k) = 120TTv{LL)k) 



On the other hand, if the rough simple pendulum model is assumed, it should be 
fuUfiled 



(6.4) 



X^w{ujk) = 1207rw(4) 
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30 40 50 

# of iterations 



80 



Figure 6.3. ur convergence. Algorithms 7 and 8. 



Besides, it is known, from fact 2, that A 
applying algorithm 3: 



27ri. So, it is decided that, when 



(6.5) 



-Mk 



(6.6) 



J^Mk 



1207r 

27r« 




-27ri 

1207r 







but a normalization constant. When applying algorithm 4, these ones are the initial 
values of Euk and Tmu- 

It is desired to compute the electromechanical mode corresponding to the gener- 
ators 1, 2,3, 8 and 10 oscillating againts the 4, 5, 6, 7 and 9 (the East side against the 
West side). As the rcducid matrix Ar + H{\) is a 10 x 10 matrix (because each Smh 
and J-"Mfc is a vector and there are 10 generators), it is needed to select the relevant 
eigenvalue and eigenvectors resulting from the factorization of the reduced matrix (see 
Algorithm 3 and 4). The chooscn procedure is as follows: 
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Q 



L 



® 



'^ & & 



Figure 6.4. New England Electric Power System 



It is defined an "objective" ao 



(6.7) 



.o = [l,l,l,-l,-l,-l,-l,l,-l,lf 



which represents the generators oscillating as described above. After performing the 
^^^Ar + H{^^^X) eigenanalysis it is obtained the matrix spectrum ■'Afe and right 
eigenvectors ^ak, where k = 1, . . . , 10 (the generator number). The products 

(6.8) pk = al °4, ^~^£m 'au 

are computed, and it is selected the mode which maximizes \\pk\\'- 



(6.9) 



arg max||pfc| 
^A = J'Afe 



'a 



oik 



'(3 = '(ik 



Note that, in the case of algorithm 3, the matrices ^£ are constant and equal to 
the initial one ^£. 



In algorithm 4, the vector ^vm is computed in each iteration by using (G.14, G.22, 



G.56 , G.57 ), and also ^vm by using analogous formulae. The matrices £ and T are 
unpdated according to: 



(6.10) '£Mk - 'vMk 

(6.11) ^^Mfe = '■'WMk 

but, possib ly, a normalization constant. 

Fi gur e 6.5 shows the modulus and phase of the right and left eigenvectors, and 
figure 6^ the modulus and phase of the 5k and uJk components of the eigenvectors. 
The ve rtic al dotted lines seggregates the variables belonging to different generators. 
Figure 6.7 shows the error evolution. The computed eigenvalue is —0.2617 + 6.4017z. 
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Figure 6.5. Right and left eigenvectors 



7. Conclusions. In this article a new approach for solving the generalized eigen- 
value problema has been introduced. The introduced algorithms can make efficient 
use of physical information regarding the shape of the sought eigenvectors. 

Appendix A. Proof of the main results.. 



From (2.1) and (2.7) it is obtained 



(A.l) XE£a + XEz = A£a + Az 

Premultiplying ( [A1| ) by T'' 

(A.2) XT'^EEa + XT^Ez = P A£a + T'^ Az 

But T'^Ez = and T'^EE = /„. So 



(A.3) 



XoL = PAEol + PAi 



Note that, but the last term T"^ Az^ this equation is an eigensystem of the n x n 
matrix !F'' A£. This matrix is usually much smaller than A. 
Introducing the matrices defined in (2.14): 



(A.4) 



P = /,„ - ESJ^^E = /,„ - Q 



These matrices are idempotent ones {V —V^,Q— Q^). Premultiplying ( |A.l| ) by 
V, and taking into account that 
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Figure 6.6. Right and left eigenvectors. 5 and ui components. 
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Figure 6.7. Error evolution. 
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(A.5) 
(A.6) 



it yields 



VEz = Ez 



(A.7) 
So, 



\Ez = VA8a + VAz 



(A. 



{\E~ A+ QA)z = VA£a 



From this equation, it is possible to solve z in function of a. In the same way, it 
is obtained 



(A.9) 



y 



\\E - A + AQ) ^ P'^ T^ AV 



Taking into account that Qz — and y^Q — Q, these formulae can be written in 
may different ways. The most symmetrical one is: 



(A.IO) 
(A.ll) 



{\E^A+[Q,A]^)z = VA£a 
y\XE -A+[Q, A].) = p^T^AV 



Therefore, from (A.IO): 



(A.12) 
(A.13) 



z = {XE-A+[A,Q]^] ^VA£a 
2/t = p^T^AV {\E-A+ [A, Q] + }"' 



assuming that the inverse matrix exists. Sufficient conditions will be provided in 
appendix P. On the other hand, it is easy to check that Vz = z, y^V = y^ . Therefore 



(A.14) 
(A.15) 



z = P{XE-A+[A,Q]^} ^VA£a 
yt ^f]^j^^AV{\E-A+[A,Q]^y^V 



and equation (A. 3) becomes 



(A.16) Xa = J'^A£a + J'''AV{XE-A+[A,Q]^} ^ VA£c 

(A.17) =Arra + H{X)a 



Note also, that by solving z from (A.8), it is obtained: 



(A.18) 



H{X) = T^AV {XE~A+ QAy^ VA£ 
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and by solving y from (A.9): 



(A.19) 



H{\) ^ Pap {\e - a + aqy^ vas 



These ones are just some few of the many equivalent ways to write H(X). Some 
of then many be more amenable for computation than (2.13). 



Appendix B. The relationship with "classical" SMA. 



The formulae (2.12.2.13) are just the generalized version of the "clasical" SMA 
(Selective Modal Analysis) formulae, as defined in section 1. 

To show the relationship, consider the problem (1.1). Let also assume that A, 
V, w are partitioned according (1.3,1.6). Therefore, the relevant subspaces in the 
generalized version are: 



(B.l) 



So 



r 1 


. 


• ^ 1 









1 . 


. 







. 


. 1 





■ /„ " 







. 


. 




_ 


. 


■ . 







T 



rr -^^rz 



^zr -^^zz 



(B.2) Pas - [ /„ ] 

as suggested by the notation. On the other hand 



4 4 



In 




j^r 



(B.3) 

Therefore 



£P ^Q^ 



In 





[ /„ ] = 



/„ 




(B.4) 

Then, 



[AQl 
' A. 
A 



rr 


Ar. 1 




r /„ 


■ 


+ 


r /„ 


■ 




zr 


A,, 




















1 






Ar. 1 






zr 





+ 













A. 
A, 



J^r 



Ar 
A., 



A., 



(B.5) 
(B.6) 






A + Arr 

A - A,, 

(\ + A„y^ 

(A-A,,)-^ 
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Besides 



(B.7) 
So 



V = Im-Q 







VA£ = 



(Bi 





U -i-m — n 





A,, 





-^rr 


Ar. 1 




\ In 1 




A.r 


A.. 









-^rr 






Azr 









(B.9) 

And, finally 



T^AV =[ln ] 

= [ Ar,\ 









Im~ 



h{x) = (t^e'p){x-a + [a,q]^} Uvas 



= [ Ar,] 



(B.IO) 





Arz (X — Azz) Azr 





[X-AzzY 



A,, 



Appendix C. Proof of theorem 

Proof. Let us recall ^e = ■'A — A. A is in the spectrum of Aj-r + H{X) and, because 
algorithm 3, ■'A is in the spectrum of Arr + H{^~^\). But 



(C.l) 



Arr + H{^-^X) = Arr + H {X) + ^^^ ^'h + o{'-\) 

oX 



By appliying a well-known perturbation formula: 



(C.2) 

Let us define 



^A- A = J< 






i-i 



e + oi^-'e) 



(C.3) 



/3t 






/3t. 



If A and ^ ^A are close enough, it is possible to neglect the higher order terms 
o(e-'^^). Then, it is clear that the algorithm converges (e-' — > 0) if and only if | p^^ \< 
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On the Other hand, from ( |2.13D , ( |A.14D and ( |A.15| ): 



(3J'^AV{XE-A+[A,Q]^} ^ E {XE - A+ [A,Q]^} ^VA£a 
(C.4) y^Ez 

Therefore 



(C.5) 



V^Ez 



But, premultiplying (2.7) and posmultiplymg (|2.9|) by Q: 



(C.6) 



(C.7) 



Qw =- Qf a + Qz 

= E£{T^E£)a + EE{T^Ez) 
^EEa 

= I3\T^E£)T^E + {v^E£)T^E 
= P'^T^E 



So 



(C.8) wtg^, = ^tgg„ = ^tjrt_£;^£:ct = /3tQ, 

On the other hand. 



So 


(CIO) 


and 


(C.ll) 


D 



w'^Ev ^ [fi^T^ + y^)E{£a + z) 

= (3^T^E£a + p^T^Ez) + (yt^£)a + y^^z 

= w^ Qv + y'^Ez 



y'lEz = wHE~ Q,)v 



LitQw 



w^{E- Q)v 



Appendix D. Proof of theorem 3.3 
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Proof. The invariance with respect to transformations £ <— £ + (/„j — E)C shall 
be proven, being the other case esentially equal. 

Firstly, note that Q and V are invariant under the transformation 

(D.l) E{£ + (/„, - E)C) Pe = EET'^E + £;(/„, - E)CP E = EET^E 

which proves the Q invariance. Ka V — I„i ~ Q, V is also invariant. Let us consider 
now: 



(D.2) 

(D.3) -{-A 

{DA)-{\E- A 

(D.5) 

(D.6) 



VA (/„ -E)C = VA (/„, -E)C 
QA} (/„, -E)C^VA {I,n -E)C 
QA] (/„ -E)C = VA (/,„ -E)C 

^{Im-E)C = {\E-A + QA}-^ VA (/„ -E)L 
~{I^-E)C=V{\E-A+ QA}-^ VA (/,„ - E) L 
-T'^A (/„ -E)C= T^AV {XE-A+ QAy^ VA (/„ - E) C 



(D.7) 

It has been used that V 



Irn - Q to go from (D.2) to (|D.3|), that E{I„^ - E) ^ 



-E) 



to go fro m (|D.3| ) t o (|D.4D and that that V{h 
go from (D.5) to (D.6). 

On the other hand under the considered transformation A„ 



{Irn - Q)ilrn - E) = En - E to 

H{X) becomes 



(D., 



(D.9) 



T^A {£ + [Ira - E)C) + T^AV {\E-A+ QAy^ VA {£ + (X,„ - E)C) = 
) T'^AE + T^AViXE- A + QAy'^VAE 

T'^A{I„, - E)C + T^AV {XE-A+ QAy^ VA{I,n - E)C = 

T'^AE + PAV {XE - a + QAy^ VAE 



It has been used the if (A) formula ( A. 18 ) to write ( D.8 ), and (D/?) to simplify 



it. 
D 



Appendix E. Proof of theorem 3.2 



The net effect of steps 3.4 and 3.5 is to impose 



(E.l) 



W^^'+^v = (^ + ^VN{'X,'QyVA) ■ 



This is because matrix ^E is proportional to the vector ^v. ^J\fy is a normalization 
constant. In fact, algorithm 4 is invariant with respect to arbitrary normalizations of 
vectors ^v and ^w (only the normalized matrices ^E, ^ J^, ^ Q and ^V are requi red). 

In order to fix the normalization ^Afv, let us focus in equations (2^) and ( A. 14 ). 



(E.2) 



v={l + ^VN{X, 3 QyVA) 'v 



u is a right eigenvector. Then, the normalization constant ^J\fy is implicitly chosen 
by enacting that this formula is valid for all i, with the same eigenvector (i.e., with 
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the same phase and absolute value). Of course, there are analogous formulae for the 
left eigenvectors. Furthermore, it is required that 

(E.3) w'fEv = 1 



Now, from (D^) it can be deduced that 



(E.4) ^MuE^+^v ={E + E^VN{^\jQyVA) E^v 

(E.5) Ev^{E + E^VN{\,^QyVA) E^v 

(E.6) ^E^v + E^z 

The last line uses the error ^ z = v — ^ z, 

(E.7) ^z ^ ^VN{XyQyVA>v 

(E.8) ^N{\yQyVA'v 



The last equation follows from ( A.12| ). Let us define the eigenvalue error ^e 
A — ^A. It is easy to check that 



(E.9) ^^^dX^^ =-N{\ ' Q)EN{\ ^ Q) 



9A2 



(E.IO) ^-^ = 2N{XyQ)EN{XyQ)EN{XyQ) 



Then substracting (E.5) from (E.4) 



^AfyE^+^v - Ev 
= E {^VN{^XyQ) - ^VN{XyQ)} E^VA^v 
^ E {^'^VN{Xy Q)EN{Xy Qy e+ 

^rN{X,^ Q)EN{Xy Q)EN{X,^ Qy e^ + . . .} E^VA'v 

^E{-'VN{xyQy€+ 

^VN{X,^Q)EN{XyQye^ + . . .} EN{Xy Q)E^VA^v 
^ E {^'^VN{X,^ Qy e + ^VN{Xy Q)EN{Xy Qy e^ + . . .} E^ 
= E {-^VNiXyaye + ^VN{X,^Q)EN{XyQye^ + ...] ^VE^ 
= E {-E^VN{XyQyVE^e+ 

E^VN{Xy Q)EN{Xy QyVEh"^ + . . .}^VEh 
(E.ll) ^ER,{^Qye)E^z 

In order to derive the last three equations it has been used that ^V^ z — ^ z and 
that ^VE = (yVE) . The last equation is the definition of matrix R^. This matrix 
is well defined so long as the matrix N{X,^Q) is bounded and ^e is small enough. 



z 
^ z 
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Sufficient conditions shall be discussed later on. On the other hand, it is obvious 
that: 



(E.12) Rz{'Q,'()^'VR,{'Q,'e) = R,{'Q,'eyV 

Let us denote by Q the matrix formed with the eigenvectors 

(E.13) Q = Evw^E 

Then 

(E.14) ER,{^ Q,^ e)Q^ z = ER,{^ Q,^ eyVQ' z 

= ER,X^Q,^eyVEvw'^E^z 

= {w^E'z)ER^['QyeyVE{^v + 'z) 

= {w^ E' z) ER,[' Qy eyv& z 

= [w'^ E^ z) ER^{^ Qy e)Eh 

So 

(E.15) ER.yQye) {E~Q)^z = [l - {lu^ & z)] ER,yQye)Eh 

On the other hand 

^^fyE^+^v -Ev= WyE {v - ^+^z) - Ev 

= {^Af,,~l)Ev-^M^E^+^z 
(E.16) =ER,yQye)E^z 



So, premultiplying by i? — Q, taking into account that (E — Q)v — and (E.15), 
it yields 



(E.17) [E-ay+'z 



,+!„_ -{E-Q)R.yQy^)iE~Qyz 



W„ [1 - (wtSJ^)] 



Premultiplying (E.ll) 



(E.18) 



^Af^w''E^+^v = 1 + w''ER,yQye)E^z 

^ , w^ER,yQye){E-Q)- 



[1 - (wt£;iz)] 
But, 

(E.19) w^E^V ^ yw^ + ^y^)E^V = ^y'^E^V 
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In this formula, it has been used the error ^y — w — ^w. The dual equations of 



(E.7,E 



(E.20) 
(E.21) 



Besides 






(E.22) 
So 



Jytgjp ^ (Jyt£;„)(u,t^jp) = {iy'^ Ev)(iy^ E'V) 



(E.23) 
and 



i„t 



[E - QyV = [1 - i'y^Ev)] h/E^V 



(E.24) 



So, as 



E^r = 



'y^E-QyV 
[1 - (iy'iEv)] 



w^'ER^i^Q^e) = w^E^VRz{^Q,^e) 

'yHE^QyVR,{^Q,^e) 







[1 - i^y^Ev)] 


(E.25) 




'yKE-Q)R.i'Q,'e) 
[1 - i^y^Ev)] 


equation 


(^.1^) 


yields: 


(E.26) 




,.r...t^,+i.._l 1 ^yHE-Q)RA^Q,^e)iE-Qyz 


But 






(E.27) 




1 - w^i;^'^ = w^Ev - w^'E^z ^ w^E'v 


and 







(E.28) 



1 — ■'y' Ev = w Ev — ■'y' Ev = ■'w Ev 



On the other hand, from (§/]) and ( |E.20D 



(E.29) 



'w'^E^z^O^^y'^Eh 
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So, 



(E.30) 



•^ w* Ev = •'vrEi^-'v + •' z) = -^ w* E-' v = w* E-'i 



Therefore, from (E.26) 



(E.31) ^Ny'^+^w^E^+^v = 1 + 



'yHE-Q)R,{^Q,^e){E-Qyz 

{iw^Eivf 



Subtituting in (E.17) 



1 



(E-32) 7Tu-^^J-T7iE - Q) 



J+i,- 



^+^wE^+^v 



Let us define the vectors 



1 I 'i/t(g-g)fl,(jg,Je)(£-Q)Jz JwE3v 



(E.33) 
(E.34) 



1 



iwE^v 



{E-Q)^z 



Jyt = 



1 



iwE'Jy 



hj^ [E - Q) 



Then, equation (E.32) can be written as 



(E.35) 



j+ 



1, -{E-Q)R.{'Q.'e) ^^ 



l + Jyti?,(JQ,J£)J5 



There is also a dual equation 



(E.36) 






N ow, le t us write the partitipation factor ■'p in terms of these new error vectors. 
From (IC.llD 



(E.37) 

Taking into account that 






(E.38) 

it is easy to check that 



'Q = 



1 



Jw'^EJv 



E'^v^w^E 



(E.39) 



'p-' 



^w^E- QYv 
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But 



(E.40) ^wt(£; - QYv = {w^ - '^y^) [E - Q) {v - ^ z) = ^ y^ {E - QY z 

because w^{E — Q) — Q and {E — Q)v — 0. Besides, 



(E.41) 
So 



iw'<Q^v = {^w'<Ev){w^E^v) = {^w^E^vf 



(E.42) 



(Jw'^E^v) 



(Remenbcr that {E — Q)^ = E — Q). Convergence results will be proved from 
formulae ( |E.35D , ( |E.36D and (^A2\) . 

Some bounds will be derived. The 2-norm, denoted as |1 • |1, will be used in the 
sequel. Firstly, note that 



(E.43) ||S-Q||<1 

In fact, let us consider an arbitrary vector v. As i? is a projection matrix 



(E.44) 



\Ev\\ < \\v\\ 



On the other hand, v can be decomposed in a component lying on the eigenvector 
V and a perpendicular component vj_: 



(E.45) 
So 



V — av + bv\ 



(E.46) ||(£;-Q)i-|| = \\a[E-Q)v + b[E-Q)vx_\\ = \\bEvi_\\ < \\bv^\\ < \\v\\ 
Let us assume that 



(E.47) 
Then 



\\'~Z\\ < &^ , VVW < Sy 



(E.48) 

As \\E-Q\\ <1, 
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(E.49) 
(E.50) 



'z|| = \\v-^v\\ <S,\\^w^E^v\\ 
\\w-^w\\ < SyW^w'fE^vW 



\\'y\\ = 



As, because of (E.29) 



(E.51) 1 = w'^Ev = Pu>^ + iy'^)E{iv + ^ z) = ^w^E^v + ^y^E'z 

it is fuUfiled 



(E.52) 



\Yw^E'v\\ < l + \\^y^E^z\\ 



Therefore, from (E.49) and (E.50), 



(E.53) 



W^w^'E^vW < l + S^SyW^w'^E^vW 



so 



(E.54) W^w^E^vW < 1 + Vl 4J,^, ^ ^ 



And (|1|) and ( |E.50|) yield 



(E.55) 
(E.56) 

So 

(E.57) 
(E.58) 



IP'zll = \\v-^v\\ <5, 
ry\\ = \\w~^w\\<5y 



\\E^v\\ < \\Ev\\ +5, 
\\^w''E\\ < \\w''E\\+Sy 



Up to now, eigenvectors v and w are arbitrary subject to the condition w^ Ev = 1. 
In the sequel, they will be chosen such that 



(E.59) 
Then, 



K = \\Ev\\ = WEwm , w''Ev = 1 



Q-^Q = Evw^E - E^v^w^'E 

= E{v- ^v) w^E + E^v {w^ - ^w^) E 
(E.60) = E^zw'^E + E^v^yE 



So, being m the dimension of Q, 
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(E.61) 



IQ-JQII = \\E^zw'<E + E^v^yE\\ 

< \\Ehw''E\\ + \\E^v^yE\\ 

< w? ^^E'^zWWiu^EW + \\EH\\\YyE\\ 



Let us consider the definition of tlie matrix i?z("'Q,-'e). It is clear, because of 
continuity with respect to ^e, that 



(E.62) 3(5,1 > 0\|P'e|| < (5,i ^ ||i?,PQ/'e)|| < \f2\\E'} V N {\' QY V E\\\y e\ 
But 



(E.63) 
So 



\\E^VN{\^QyVE\\ < \\E^V\\\\NiX,^Q)\\\\^VE\\ 

<\\E-=Q\\\\NiX,^QmE-^Q\\ 
< \\NiX,'Q)\\ 



(EM) 



3Sei > 0\|P'e|| < S,i ^ ||i?,PQ,^e)|| < V2||7V(A/Q)|||P'e|| 



Because of (E.61), Q and ^ Q are close whenever ^ z and ^y are small. Therefore, 



(E.65) 36,1 > OMPzIl < ,5,1 A |Py|| < ,5,i ^ ||A^(A,^Q)|| < V2||iV(A, Q)\\ 
So, 



(E.66) 



35ei > A 36,1 > OWe\\ < 6,i A \y ~z\\ < 6,i A |Py|| < J.i 

||i?.(-'Q,-''e)||<2||iV(A,Q)|||P'6|| 



Bounds on N{X, Q) will bw provided in the next appendix. Besides, from ( E.35| ) 
and (|E.36D 



(E.67) r+'~z\\ < [l + 2\\^y^R,{^Q,^eyz\\]\\E-Q\\\\R4^Q,^e)\\\\^S\\ 

(E.68) |P"+iy|| < [l + 2|Pyti?.eQ,^eF^II]l|i?- QIIP.es, ^■e)IIIPy|| 

whenever |py f i?2(^Q,^e)-'z|| < 1. Therefore, 



(E.69) 



36,1 > A 3,5,1 > 0\|P'e|| < ,5,i A \\' S\\ < 6,i A \\'y\\ < 6, 



IP'+i^ll < 2 
IP'+^y|| < 2 



l + 4||7V(A,Q)||,5,i52^ 
l + 4||iV(A,Q)P,i,52i 



|A^(A,Q)||||-''e|||P'5|| 



GENERALIZED SELECTIVE MODAL ANALYSIS 



33 



Let us define 



(E.70) 
Then, 



<522 = inm{6zi, 



A\\N{X,Q)\\d, 



,}>0 



(E.71) 3S,i > A 36,2 > 0\\\^e\\ < d,i A \\^ S\\ < 5,2 A |Py|| < 6,2 

\V+'~z\\ < 4||iV(A,Q)||||J-6|||p-z~|| 
r+'y\\ < 4||iV(A,Q)|||K6|||Ky|| 



On the other hand, from ( |C.2[) , it is clear that 



(E.72) 36,2 > 0\|P"e|| < 6,2 => |P"+^e|l < 2\\ = p-'We\\ = 2\\^SWyW^ 



because of ( |E.42| ). Let us define 



(E.73) 



6,3 = inm{6ei,6,2} > 



Then, from (|E.71D and (|E.73D , 



(E.74) 36,3 > A 36,2 > 0\|Pe|| < 6,3 A \\^ S\\ < 6,2 A \\^y\\ < 6,2 

IP+i^ll < 4||iV(A,Q)|l|N|p5|| 
r+'y\\ < 4||iV(A,Q)|l|N||py|| 

Ip'+i^ii < 2\rz\\\\'y\W4 



Let us now to define 



(E.75) 



(5,4 = mm{6,3, giijvfA Q)\\ ^ ^ ^ ' '^^^ "" ™n{(5^2, ^l 



Then 



(E.76) 



36,4 > A 3,5,3 > OMPell < ,5,4 A \\^ S\\ < 6,3 A \\^y\\ < 6,4 

\'^'y\\ < Wy\\ 



which imphes that the algorithm converges if inizializated in a neighborhood of the 
solution {X,v,w). In fact, if pzH < pu''''i?-'w||^23|| the n, bec ause of ( E.33 ) |pz|| < ,5,3. 
There is an analogous formula for |py||. In that case, E.76 applies and the algorithm 
converges. 



In order to asses the quotient convergence factor, let us write E.74 as 
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(E.77) 



log|p'+i5 
log||J"+iy 
log|p+ie 



< 



1 1 

1 1 

1 1 1 



logll-J^II 

logll^yl 
loglpejl 



log4||7V(A,Q)|| 

log4||7V(A,Q)|| 

log 2 



By substituting < by = a majorant sucession is obtained. Its asymptotic be- 
haviour is controlled by the dominat eigenvalue of the state matrix, which happens 
to be 1 + \/2- Therefore, it is expected that, asymptotically 



(E.78) 



log|p■+le|l<(l + ^/2)log|P"6 



which completes the proof of the theorem. 

Appendix F. Some bounds on N{\, Q). 

The purporse of this appendix is to provide some bounds on the matrix iV(A, Q), 
when A is an eigenvector of the pair [E, A) and the matrix Q is formed from its right 
and left eigenvectors: Q — Evw'^E. 

It shall be assumed that the pair {E, A) is solvable 0, that is, the pencil ^E — A is 
regular for all /^ but a finite number. Of course one of these /i is the sought eigenvalue 
A. Moreover, it shall be assumed that A is a single eigenvalue. 

Let us assume that 



(F.l) N{\Q)-^ ^\E-A + AQ 

is singular. Then, there is a vector x such that 

(F.2) {\E^A + AQ)x = Q 

However, it is possible to write 

(F.3) x = av + z ,w'^Ez^Q ,w^Ev^l 

So Qz = 0, Qu = Ev. Then, 

(F.4) aXEv + {\E~ A)z^Q 

Premultiplying by the left eigenvector w, 



(F.5) 



aA = 



Then, from (F.4), it is obtained that {XE — A) z — 0. But this is impossible, 
becuase it is assumed that A is a single eigenvalue. Therefore, the matrix \E — A-\- AQ 
is regular and A^(A, Q) is bounded. 

A similar reasoning can be done if 



(F.6) 



N{\,Q)-^ ^XE~A+QA 
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(F.7) Ni\,Q)-' =XE-A + AQ+QA 

Then, it is obtained that 

(F.8) a (XE + QA)v + {XE - A+ QA) z ^ 

But Av — XEv, so 

(F.9) 2aXEv + {XE- A+ QA) z = 

Premultiplying by w^ E, and taking into account that w^ EQ = w^ E, it is obtained 
that 



(F.IO) 

So, from (FJ), 



2aX + Xw^'Ez = aX = 



(F.U) 



{XE-A+QA)z = 



which implies that the matrix XE — A+ QA is singular, or that z = 0. As it has been 
proved that XE — A + QA is regular, it must be z = 0. Then, the singular vector 
X — V. From (|F.S| ), it must be XEv = 0. So: 

• If A = 0, then N{X, Q) is singular, and i; is a singular vector. 

• If A 7^ 0, then we must have Ev = 0. But Ev = jAv ^ 0, because if Av — 0, 
then A = 0. So, the matrix N{X, Q) is regular. 

Appendix G. Composite models. 

Let us define: 



(G.l) 



Xm 
XI 

xo 

XA 



The system equations are 



(G.2) 







' di&g{Ei,E2,. 


■.El) 





" 




XM 



















XI 



















Xo 














_ 




XA 


diag(Ai,^2,- 


■ ,A,) 


diag(Bi,B2,.. 


■,Bi) 










XM 


diag(Ai,^2,- 


■ ,Ai) 


diag(Bi,B2,.. 


■:Bl) 


/ 







XI 







-Jii 




/ 


-J12 




Xo 







-J21 







^J22 




XA 
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where / denotes the identity matrix. 

This hnear dynamic system has eigenvectors which can be partitioned analogously 
to the variables and equations. Specifically, it is fulfilled: 



(G.3) 

(G.4) 
(G.5) 
(G.6) 
(G.7) 
(G.8) 

(G.9) 
(G.IO) 



XEvm — Avm + Bvi 

— Cvm + Dvi — vo 

= -Jiwi + Wo - J12VA 

= -J21V1 - J22VA 
Xw\,jE = wIjA + w\C 

— w\jB + w\D — (wqJh + 110^321 



= -wUJu 



w\j22 



By introducing the vector 



(G.ll) Wq = WqJii + ■w\J2i 

and writing the previous equations in a subsystem basis, it is obtained: 



(G.12) 
(G.13) 

(G.14) 

(G.15) 
(G.16) 

(G.17) 



XEkVMk — AkVMk + BkVjk 
vok = CkVMk + DkVik 





Jii 




J12 




VI 




Vo 






J21 J22 




. '"A 









^wli^Ek = wlj^Ak + wlj^Bk 


^hk = ^lik'^k + w\^Dk 


I 


ulw\ 




■ Jn 




J12 ' 
J 22 


= 


w^Q 





It is easy to check that 



(G.18) 



T'ES = /,„ =^ T.f.Ek£_ 



-Mk 



I, 



■mk 



\Jk 



where mk is the number of columns of Euk- Also, the formulae (pJ-2.10) can be 
written as: 



(G.19) 
(G.20) 
(G.21) 

(G.22) 
(G.23) 
(G.24) 



VMk — ^MkCtMk + ZMk 
Vik = Zik 
VOk — ZOk 
VA == ZA 



^Mk^kZMk — 



WMk = ^Mk^Mk + VMk 
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(G.25) 
(G.26) 
(G.27) 
(G.28) 

and 



wik = yik 
wok = yok 

WA = VA 

8l,j^.EkyMk = 



(G.29) a 

From (|G.12|) and (|G.19[]G.2l[ ) 



OiMl 




/3m 1 


aM2 


/3 = 


I3m2 


CUMl 




L PmI 



(G.30) XEk^MkOiMk + XEkZMk — AkEMkO-Mk + A^ZMk + BkZjk 



Premultiplying by ^j^jf., and taking into account ( G.23 ) 



(G.31) Aajv/fe — ^Mk^k^MkO^Mk + ^Mk^kZMk + ^Mk^k^Ik 

Let us define the projections Vk y Qfc by: 

(G.32) Vk — Imk — EkEukJ^Mk^k = Irak — Qk 

It is easy to check: 



(G.33) 
(G.34) 



VkEk^Mk — 

VkEkZMk — EkZMk 



Taking the above equations into account, and premuhiplying (|G.30 ) by Vk 



(G.35) 



So 



XEkZMk = VkAkEMkOLMk + VkAkZMk + VkBkZjk 



(G.36) 



{XEk - Ak + QkAk) ZMk = VkAkEMkCtMk + VkBkZik 



Taking into account that VkZuk = ZMk and QkZuk = 0, the above equation 
yields: 



(G.37) ZMk^Vk{\Ek-Ak + [Qk,Ak\ + ) ^Vk{AkEMkaMk + BkZjk} 



After substitution in equation (|G.3l|) it yields: 
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(G.38) XaMk = (Ark + HAkW) auk + {Brk + HekW) zik 

where the following matrices have been defined: 



(G.39) 
(G.40) 
(G.41) 
(G.42) 



A-rk — J'MkAk^Mk 



Br 



J^Ml-Bh 



HAkW = ^lik^kVk {XEk - Ak + [Qk,Ak]+y' VkAk£_ 
HBkW - ^likAkVk {XEk ~Ak + [Qk,Ak] + r^ VkBk 



Mk 



On the other hand, from equation (4.11), it is obtained 



(G.43) 



ZOk — Ck^MkCtMk + CkZMk + DkZjk 



After substitution of (G.37) it is obtained 



(G.44) zok = [Crk + Hck) aMk + [Dk + Hok) zik 

where the following matrices have been defined: 

(G.45) Crk = Ck^Mk 

(G.46) Hck{X) = C'kVk {XEk - Ak + [Qk,Ak] + r' VkAkSuk 

(G.47) HokiX) = CkVk (XEk ^ Ak + [Qk,Ak]+r' VkBk 

Let us define the matrices: 



(G.48) 
(G.49) 
(G.50) 
(G.51) 
(G.52) 
(G.53) 
(G.54) 
(G.55) 

It is easy to check that Ar — T^AE, as computed according the general formula. 
Therefore, the proposed notation is consistent. Note also that all these matrices are 
diagonal-block matrices, which eases its computation. Then, equations ( G.38 , p~44 ) 
can be written as: 



Ar 


^ diag(^ri • 


■Arl) 


Br 


== dl&g{Brl . 


■Brl) 


Cr 


= diag(Cri . 


■Crl) 


D 


= diag(Di . . 


• A) 


Ha 


= Aia.g{HAi ■ 


■ ■Hai) 


Hb 


= dia.g{HBi ■ 


■ ■Hbi) 


He 


= diag(iJci ■ 


■ ■ Hci) 


Hd 


= diag(iJ_Di . 


■ ■Hdi) 



(G.56) 
(G.57) 



Xa = [Ar +HA{X))a + {Br + Hb{X)) zj 
zo = {Cr + Hc{X)) a + {D + Hd{X)) zi 
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Xa = {A,. + Ha{X)+ 
(G.58) [{Br + Hb{X)) 0] 

Therefore 



Jn-{D + HDiX}) Ji2 
J21 J22 



Cr+HciX) 





H{X) = Ha(X) + [{Br + Hb{X)) 0] 
(G.59) 



Jn-{D + HD{X)) J12 
J21 J22 



Appendix H. Proof of theorem 5.1 

Let us assume that 



Cr+Hc{X) 





(H.l) 
But 



[A - J-^A {E - Q)] ^v = ^XQ^v 



(H.2) 
So 



Q^v = {E£T^E){8'a + ^ z) = E£'^ 



a 



(H.3) [A - ^-^X {E - Q)] (£^a + ^'z) = ^XES^a 

Premultiplying by J-"^ , and taking into account 



(H.4) 
it yields 



^^ iE-Q)= T^E - PeEET^E = 



(H.5) T'^APa + T^A^z^iX^a 

On the other hand, as 



(H.6) 



{E^Q)E^EE- EEPeEE = 



equation (H.3) yields 



(H.7) 



A£^a + [A - ^'-^A {E - Q)] h. = ^ XE£h 



Premultiplying by V , and taking into account that V{E — Q) = E — Q, that 
VE£ ^{E- Q)£ = 0, and that Q^ z = 0, 
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(H.; 



VAE^a + VA^z = ^-^\E^>. 



This is, esentially, formula (A. 7). So, it is possible to conclude (A. 12) 



(H.9) ^z^V{^-^\E-A+[A,Q]^] ^VAE^a 

and, by substituting in (|H.5|), 



(H.IO) 



^\^a= {Arr + H{^-^\))^a 



Analogous reasoning can be done by using the left eigenvector, which proves the 
theorem. 



Appendix I. Proof of form ula ( 5.26 ). 

From the ^V definition ( 5.17 ): 



(1.1) ^V =[A-^-^\{E-EEPe)] E£ 
By using the Shcrmann- Morrison lemma: 

^V = [{A - ^-^XE) - ^~^\E£pE] "^ E£ 

= Ua- 3-^xEy^ - J-^A {A - ^-^xEy^ 

In + T^E {A - ^-^XEY^ E£^-^X\ ^ Pe {a - ^-^XE) "^ \ E£ 

(1.2) = \^V - '-^X^V [in+'-^XT^E^VY^ T^E^V^ 
Taking into account that W^^ = T'^ EW and that ^M = T'^E^V, 

W"^ = Pe^v 

= ^M - ^-^X^M [In + ^^^X^M] "^ ^M 

= ^M [In + ^-^X^M] "^ ([/„ + ^-^X^M] - ^-^X^M) 

(1.3) =^M[ln + ^'^X^M]~^ 
So, 



(1.4) 



W = [In+^^^X^M] ^M-^ = ^M-^+^-^XIn 



Appendix J. Proof of theorem 

Let us consider the equation 



(J.l) 



[A - X{E -Q)]v = E£a' 
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(J.2) 



V = £a + z , T* Ez = 



Premuhiplying (J.l) by J^^^ 



(J.3) T^A-X{E-Q)]v = a' 

As J^^(S-S) = 0, 



(J.4) 



T'^ASa + T^A 



z — a 



On the other hand, from (|j.l[), as {E — Q)£ — 



(J.5) 



Av - \{E - Q)z ^ EEa' 



So, 



(J.6) A8a +[A- \{E -Q)]z^ E£a' 

AsV{E- Q) = E- Q, VE£ = and Qz = 0, premultiplying by V, 

(J. 7) VA£a + VAz-\Ez^O 

And, as Qz = and Vz — z, 



(J.8) 



z = V[XE-A+[A, Q]+]~^ VAEa 



So, from (J.5), 



ol = pAEa + PAV \\E-A+ \A., 0]+]"^ VAEa 
(J.9) = (A^r + iJ(A))a 

On the other hand, let us define 



(J.IO) 
So, 



V = \A-\{E-(2)Y^EE 



(J. 11) v = VoL =Ea^z 

Premultiplying by J-^E, and taking into account J-^EE — In, 



(J.12) 



J^^EVa' = a 
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So 



(J. 13) {T^EV) a^Af: 



-1 

a ^ a' 



Now, a' is arbitrary. So, as ( |J.9| ) and ( J. 13 ) hold for any a', the matrices must 
be equal: 



(J. 14) Af^Arr + H{X) 

which proves the claim. 
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